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Forecasts of busy season trunk group traffic loads are required for 
planning the Bell System's message network. Forecasting algorithms 
currently in use obtain estimates of future loads by multiplying the 
most recent measurement of busy season load by an aggregate growth 
factor. Because of statistical errors in measured loads and differences 
between individual trunk group and aggregate growth factors, the 
resulting forecasts can have large statistical errors. In this paper we 
extend earlier work to develop a new algorithm, called the sequential 
protection algorithm (spa), based on a linear two-state Kalman filter, 
together with logic for detecting and responding to unusually large 
measurement errors or changes in trend. In typical applications of 
Kalman filtering, the statistics of system noises, measurement errors, 
and initial conditions are known and the filter parameters (Kalman 
gains) are selected accordingly. For our application, however, these 
statistics cannot be determined without error. Consequently, we de- 
velop a method for selecting robust filter parameters which provide 
improved performance, independent of system noises, measurement 
errors, and initial conditions. In particular, under the assumption of 
linear growth for 5-year intervals, the average rms 1-year forecast 
error of spa is about 10 percent less than that of the existing algo- 
rithms. Field test results confirm the theoretical results presented 
here. Accordingly, specifications have been written for inclusion of 
spa in the Bell System 's standard trunk forecasting systems. 

I. INTRODUCTION AND SUMMARY 

Forecasts of busy season (yearly peak) trunk group traffic-loads are 
required for the planning of the Bell System's message network. These 
forecasts are used to design traffic networks which minimize the cost 
of the trunks required to satisfy anticipated customer demands. 
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The standard load forecasting algorithms currently in use in the Bell 
System obtain estimates of future loads by multiplying the most recent 
measurement of busy season load by an aggregate growth factor; for 
example, the average of the growth factors obtained by trending the 
total office loads at each end of the trunk group. Descriptions and 
comparisons of the various algorithms currently in use are given in 

Ref. 1. 

As explained in Ref. 1, these algorithms have two significant sources 
of error: (i) Because of the finite amount of data upon which measure- 
ments are based, measured loads can have large statistical errors; 
standard deviations fall in the range of about 5 to 40 percent depending 
upon load size and type of measurement system. 2 (ii) Individual trunk 
group growth factors can differ from the aggregate growth factor; 
standard deviations of 6 percent have been observed. These forecast 
errors are significant since they lead to an increase in the reserve trunk 
capacity required to satisfy customer demands. 3 

To reduce forecast error and, hence, reserve trunking capacity, a 
new algorithm, called the sequential projection algorithm (spa), has 
been developed to forecast busy season traffic loads within the Bell 
System. The spa is based on a linear two-state Kalman filter model, 
whose use in traffic forecasting was studied first by David and Pack, 1 
together with logic for detecting and responding to outlier measure- 
ments, i.e., unusually large measurement errors or changes in trend. 

As discussed in Ref. 1, David and Pack tested several Kalman filter 
models— some with as many as eight state variables and four data 
variables. In summary, for the planning interval of interest (1 to 5 
years ahead) none performed consistently or significantly better than 
the relatively simple two-state (traffic load and incremental growth), 
one-data variable (measured load) model. 

The two-state Kalman filter establishes a linear trend for individual 
traffic loads as follows: As illustrated in Fig. 1, the level, or smoothed 
base load, is a weighted average of the most recent measurement, or 
base load, and the previous 1-year forecast. Similarly, the smoothed 
growth increment is a weighted average of the measured and previously 
forecasted increments. (The measured increment is, by definition, the 
measured load minus the previous year's smoothed base load.) 

The performance of the Kalman filter, as measured by mean square 
forecast error, depends upon the filter parameters, i.e., the gains a n 
and /?„ in Fig. 1, the standard deviation of load measurement and 
growth estimation errors, and upon the assumed evolution of the true 

load. 

Since spa will be used under a variety of possible operating condi- 
tions, the filter parameters could be tuned to provide optimal perform- 
ance, i.e., minimum mean square forecast error, for each application. 
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SMOOTHED LOAD = FORECASTED LOAD (1-a„) + MEASURED LOAD (a„) 
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Fig. 1 — Basic operation of spa. 



For example, in the Bell System we employ two types of load mea- 
surement systems; Bell Operating Companies obtain load estimates 
from a direct measurement of trunk group usage, while Long Lines 
derives estimates from point-to-point, e.g., end-office to end-office data 
provided by the Centralized Message Data System (cmds). 2 For trunk 
group data, the standard deviation of measurement error is in the 
range of about 5 to 10 percent, depending on load size. For point-to- 
point data, the range is about 10 to 40 percent. (See Fig. 2.) Accord- 
ingly, different parameters could be used for different measurement 
systems and for different load ranges. 
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Fig. 2 — Sampling error vs. offered load. 
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Indeed, the initial development of spa 1 was intended for use with 
trunk group data only, and the objective of the current study was to 
extend the original work to applications using point-to-point data. As 
suggested above, one possible solution would be to develop multiple 

versions of spa. 

Instead, in this paper, we develop a single, robust spa whose param- 
eters are selected to provide improved performance over the entire 
range of operating conditions, including the use of either trunk group 
or point-to-point data. 

The use of robust parameters is important for two reasons: (i) A 
single spa for all applications should be simpler to implement and 
maintain than multiple versions, (ii) The actual values of the statistical 
parameters for each application cannot be determined without error, 
and our results show that erroneously assumed values can lead to a 
performance substantially worse than that of conventional projection 
methods. Of course, as with any robust technique, we pay a premium 
by receiving less than theoretically optimal performance for protection 
against the possibility of a performance worse than conventional 
methods. 

Section II of this paper provides a qualitative overview of the 
functions performed by spa that include procedures for detecting and 
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Fig. 3— Sequential projection algorithm. 
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responding to outliers, as well as the smoothing and projection func- 
tions of the Kalman filter. Section III defines the mathematical model 
and procedures for selecting robust filter parameters; Section IV gives 
numerical results; Section V develops the outlier detection thresholds; 
and Section VI gives the conclusions. 

II. OVERVIEW OF SPA 

As shown in Fig. 3, spa is composed of five major components: 
algorithm initialization, base-load adjustments, outlier detection, 
smoothing, and projection. In the following paragraphs, we provide a 
qualitative description of the operations performed by each of these 
components. 

2.1 Initialization 

As indicated in Fig. 4, the smoothed base load in the first year of 
operation, and in certain other cases discussed in Section 2.3, is equal 
to the measured base load. The smoothed growth increment is calcu- 
lated by multiplying this base load by a growth factor obtained, for 
example, by trending the total office loads at each end of the trunk 
group. 1 

2.2 Base load adjustments 

In the second and subsequent years of operation, spa updates both 
the level and incremental growth by comparing the most recent base 
load with the previous 1-year forecast of that same load. Since the 
base and forecasted loads must correspond to the same traffic routings, 
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Fig. 4 — Initialization of spa. 
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differences between the previous and current routings are accounted 
for by adding an adjustment term to the base load so that the adjusted 
base load agrees with previous routing assumptions. After performing 
the outlier detection and smoothing functions described below, the 
routing adjustment term is subtracted from the smoothed base load so 
that it agrees with the current routing assumptions. 

Furthermore, the previously forecasted load is adjusted to remove 
the impact of deterministic events, such as a proposed tariff change, 
that were predicted but did not occur. 

2.3 Outlier detection 

Under the assumption that the observed forecast error (i.e., the 
difference between the forecasted and measured loads) has a Gaussian 
distribution with zero mean, the linear Kalman filter upon which spa 
is based provides the minimum attainable mean squared forecast 
error. 4 In practice, however, the normal statistical errors (because of 
the finite measurement interval, day-to-day load variation, random 
variations in cmds point-to-point sample size, 2 and growth errors) are 
occasionally contaminated by wiring errors, recording errors, or unex- 
pected changes in the trend of the true load. In such cases, when the 
observed forecast error deviates from a Gaussian distribution, the 
linear Kalman filter model can have a mean squared forecast error 
which is substantially greater than the minimum. 4 

Our approach to this problem, which is based in part on the nonlinear 
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Fig. 5— Response of outliers (sign not repeated). 
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Kalman filter model described in Ref. 4, is the following: If the 
difference between the possibly adjusted base and forecasted loads 
exceeds present thresholds, the base load is declared to be an outlier. 
In response, spa will adjust the base load or restart at the measured 
level, depending only upon whether an outlier of the same sign oc- 
curred in the previous year. 

If an outlier of the same sign did not occur in the previous year, spa 
replaces the base load by the nearest threshold value as indicated by 
the vertical arrows in Fig. 5. Qualitatively, the underlying assumption 
here is that in most cases such an outlier signals an invalid or atypical 
measurement caused by, for example, a recording error and not a 
change in trend. Formally, this modification of the Kalman filter is 
equivalent to that proposed by Masreliez and Martin in Ref. 4. 

Alternatively, as indicated in Fig. 6, if an outlier of the same sign 
occurs in two consecutive years, spa restarts at the measured level. 
The assumption here is that in most cases two consecutive outliers of 
the same sign signal a change in trend. In theory, additional improve- 
ment could be obtained by restarting at the previous outlier and then 
smoothing with the current measurement. However, such a procedure 
would be more complicated to implement, and our studies show that 
it would have negligible impact on performance. 

As indicated by Huber's studies, 5 and as supported by our field test 
results, 6 adequate protection against outliers is obtained when the 
thresholds are set anywhere in the range of about one or two times the 
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Fig. 6 — Response to outliers (sign repeated). 
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rms observed forecast error. This result is important: we show that it 
allows us to use thresholds which are independent of the number of 
data points processed and the type of measurement system. Numerical 
values for these thresholds are provided in Section V. 

2.4 Smoothing 

The possibly adjusted base and forecasted loads are combined to 
produce a smoothed base load and a smoothed growth increment. As 
indicated in Fig. 1, the smoothed base load is a weighted average of 
the base and 1-year forecasted loads. Similarly, the smoothed growth 
increment is a weighted average of the measured and forecasted growth 
increments. Procedures for selecting the appropriate gains a n and /?„ 
are described in Section IV. 

2.5 Projection 

As indicated in Fig. 1, the smoothed base load and growth increment 
are combined to establish a linear projection. That is, the projected 
load for the yfeth future year is obtained by adding k smoothed growth 
increments to the smoothed base load. In practice, the resulting trunk 
group load forecasts can be adjusted to include user supplied estimates 
of the impact of deterministic events (caused by, for example, proposed 
tariff changes) or to agree with other aggregate load forecasts. 

III. MATHEMATICAL MODEL 
3.1 General 

Although spa is based on a linear two-state, i.e., trunk group load 
and incremental growth, Kalman filter, we considered more complex 
models with additional state variables, e.g., aggregations of other trunk 
group loads and growth factors. Therefore, for completeness, we first 
summarize the equations which define the general discrete-time linear 
Kalman filter. For a more complete discussion, see Ref. 7. 

The true time-behavior of the state variables is assumed to be 
defined by the linear transition equation 

X n+1 = 4>Xn+Wn+ Un, (D 

where X n is an s- vector of true state variables in period (year) n, <j> is 
an s X s transition matrix which may depend on n, W n is an s-vector 
of zero-mean random modelling errors, and Un is an s-vector of 
deterministic changes in state. 
The one-period-ahead projection formula is 

X n +l,n — <i>Xn,n + Un, (2) 

where X n+k ^ denotes an estimate of X n +k based on data through period 
n. The "smoothed" estimate X n ^ is calculated by 
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X n ,n — Xn,n-\ + K n (y „ — HX n ,n—\)> (3) 

where K n is a d x s "Kalman gain" matrix andy n is a d- vector of 
measurements in period n related to X n by 

y n = HX n + V n , (4) 

where H is a d X s matrix and V„ is a d- vector of zero-mean measure- 
ment errors. 

The optimal gain matrix K„ is determined recursively by the equa- 
tions 

K n = PnH T {HP n H T + i?)" 1 (5) 

S n = (I - KnH)P n , (6) 

and 

P n+1 = <t>S n <t> T + Q, (7) 

where R is the covariance matrix of measurement errors, i.e., R = 
E(V n Vn), Q is the covariance matrix of modelling errors, and So is an 
estimate of the covariance matrix of the initial state estimate Xo.o- 
Furthermore, it is assumed that V, and Wj have zero mean and 
are pairwise uncorrected among themselves and with each other for 
alii,/ 

As discussed in Section I, we will be concerned with filter perform- 
ance for nonoptimal gains. In this case, the covariance matrix P n +i of 
the projection error (X„ +i , n — X n +i) is related to the covariance matrix 
S n of X n ,n by eq. (7), but S n is related to P„ by 

S n = (/ - K n H)Pn(I - K n H) T + K n RK T n , (8) 

which reduces to eq. (6) only when K„ is given by eq. (5). 

3.2 Two-state model 

As noted in Section I, we considered models with as many as s = 8 
state variables and d = 4 measurement variables. However, our studies 1 
showed that none performed consistently or significantly better than 
the simple two-state, one-data variable model defined by the equations 
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and 

y n = x n + v n , (10) 

where x n and x n are, respectively, the true load and true incremental 
growth in year n, and y n is the measured load in year n. These 
equations correspond to eqs. (1) and (4) with 
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<*> = 



1 1 
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(11) 



and 



#=[1,0]. 



(12) 



The more complex models were rejected since the additional vari- 
ables that we considered, e.g., aggregations of other trunk group loads 
and growth factors, are usually not highly correlated with trunk group 
load. Moreover, when the correlation is high the improvement in 
performance is minimal and, more importantly, when high correlation 
is incorrectly assumed, the penalties outweigh the expected benefits. 1 

To complete the description of our two-state model, note that the 
smoothing eq. (3) becomes 



■ikn.n — 



Xn,n 
Xn,n 



X»Mi-l + Otniyn - X n , n -l) 
Xn*-\ + Pniyn - Xnjt-i) 



(13) 



where a„ and fi n are the Kalman gains in year n. The 1-year-ahead 
projection formula is 

(14) 



in+l,n — 



Xn,n 
X n ,n 



Note: Our analysis of the two-state model will ignore the possibility of 
deterministic changes in state; that is, we assume U n = 0. Also, by 
combining eqs. (13) and (14) it follows that eq. (13) may be written in 
the form 



Xnji 
X n ,n 



(1 - a n )Xnj,~i + a„y n 

(1 - pn)X n . n -l + (Snhn ~ X n -l,n-l) 



(15) 



which emphasizes that x n ,„ is a weighted average of the forecasted and 
measured levels and, similarly, x„, n is a weighted average of the 
forecasted and measured increments. 

Finally, note that R = a 2 , where a is the standard deviation of load 
measurement error; the dependence of a on load size and type of 
measurement system is discussed in Section 4.2.2. In the following 
section, we define the spa initialization procedure and in Section 3.4 
we describe our procedure for selecting values for the filter gains a n 
and fin and the associated assumptions about Q. 

3.3 Initialization 

As explained in Section 2.1, the smoothed-base load and growth 
increment in the first year of operation are given by 

*o,o = yo (16) 
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and 

xojo = gyo, (17) 

where g is an aggregate growth factor. For example, if g = 0.1, the 
increment is 10 percent of the base load. 

The normalized covariance matrix So/o 2 = (Sij/o 2 ) of the initial state 
estimate is given by 

s,i(0) E(xo-yo) 2 



= 1, (18) 



s 22 (0) E(gx -iyo)' 



_ E[(g-g)x + g(x l) -y )Y 
o 2 

= ^ + i 2 , (19) 

o 

where g is the true growth factor and o e is the standard deviation of 
the estimate g, and 

si 2 (0) _ E[(xo - yo)(gx - gy )] 
a a 

_ E{(xo - y )[{g - g)xo + g(x - y )]} 
a 2 

= g. (20) 

3.4 Filter parameters 

As discussed in Section I, our objective is to select robust values for 
the gains a„ and /?„; that is, values which will provide improved 
performance over the range of operating conditions. Our approach to 
this problem starts with the following idealistic assumption, which will 
be relaxed in a later section: 

We first assume that the true load displays constant incremental 
growth. That is, we first assume that Q, the covariance matrix of 
modelling errors, is identically zero. Under this assumption, it follows 
from eqs. (5) to (7) that the optimal gains depend only on the ratio 
So/o 2 . In turn, for a given value of g, it follows from eqs. (18) to (20) 
that So/o 2 is determined by the ratio 

G 2 - 7-7^72. (21) 

(o/xoY 

In Section 4.2, we display the optimal gains and corresponding 
performance (as measured by mean square forecast error) for various 
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known values of G. More importantly, we show that an erroneously 
assumed value for G can lead to a performance worse than that of the 
conventional projection method; equivalently, since the conventional 
method is used to initialize spa, the mean square forecast error can 
increase with the number of data points processed. 

Next, in Section 4.2.3, we show that the gains corresponding to one 
particular value of G provide a performance that is nearly independent 
of the actual value of G. We call these gains the "robust gains for linear 
growth." 

Finally, in Section 4.3, we relax the assumption of linear growth and 
show that certain constant gains, derived from the "robust gains for 
linear growth," provide improved performance in the presence of 
system noise, as well as the ideal case where the trend remains 
constant. 

IV. NUMERICAL RESULTS 

4.1 Examples 

To help explain our results, we first consider several examples which 
show how the gains and the mean-square forecast error depend upon 
the ratio G. In these examples, we assume that the true load has 
constant incremental growth, and we assume two-point distributions 
for measurement and growth estimation errors. That is, the measured 
load is either high or low by equal amounts and with equal probability, 
and the initial estimate of incremental growth is either high or low by 
equal amounts and with equal probability. Moreover, in each example 
we assume that the aggregate growth factor g is zero. 
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Fig. 7— Optimal filter: no measurement error (G = oo). 
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4.1.1 No measurement error (G = *>) 

In the example illustrated in Fig. 7, we assume no measurement 
error, i.e., a 2 = 0, but an error of plus or minus one unit in the initial 
increment so that G = oo. Thus, with this information, the true load 
lies along the dashed line marked A with probability 0.5 or along the 
dashed line B with probability 0.5. Accordingly, the initial mean square 
forecast error is 0.5(1) 2 + 0.5(-l) 2 = 1. 

As shown in Fig. 7, the possible measurements in year 1 are Mi or 
Mi with equal probability. However, since there is no measurement 
error and since two points determine a straight line, it is clear that in 
either case the forecast error can be reduced to zero after processing 
the second measurement; examination of eq. (13) shows that the 
appropriate gains are <x\ = /?i = 1. 

4. 1.2 No growth error (G « 0) 

At the other extreme, suppose there is no error in the initial growth 
estimate, i.e., o\ = 0, but an error of plus or minus one unit in the 
measured load so that G = 0. As shown in Fig. 8, the trend line is 
either A or B with equal probability; hence, the initial mean square 
forecast error is 1. The possible measurements in year 1 are Mi with 
probability 0.25 [since A and a high measurement occur with proba- 
bility (0.5) (0.5)], Mi with probability 0.5, and M 3 with probability 0.25. 

Since Mi and M3 correspond uniquely to A and B, respectively, the 
forecast error can be reduced to zero after processing either of these 
measurements; eq. (13) and Fig. 8 show that the appropriate gains are 
ai = 0.5 and fii = 0. However, if M7. occurs, no new information is 
obtained. Consequently, after processing the second data point, the 



A MEASURED LOAD 
» FORECASTED LOAD 
— TRUE LOAD 

A«l 






AM 3 



1 2 

YEAR 



Fig. 8 — Optimal filter: no growth error (G = 0). 
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Fig. 9 — Optimal filter: equal measurement and growth errors (G = 1). 

mean square forecast error is 0.5(0) 2 + 0.5(1) 2 = 0.5. Thus, the rms 
forecast error is reduced by a factor of 1/V2. 

4.1.3 Equal measurement and growth errors (G = 1) 

For values of G between the above extremes, i.e., < G < oo, the 
reductions in forecast error might be expected to fall between those 
corresponding to the extreme values of G. But, this is not generally 
true. 

For example, consider the case illustrated in Fig. 9 in which G = 1.0. 
We assume uncorrelated errors of plus or minus one unit in the 
measured load and in the initial increment. Thus, as shown in Fig. 9, 
the possible trend lines are A, B, C, or D, each with equal probability. 
Consequently, the initial mean square forecast error is 0.25(2) 2 + 
0.5(0) 2 + 0.25(-2) 2 = 2. 

In year 1 the possible measurements are Mi with probability 
[0.25(0.5) = 0.125], M 2 with probability 0.375 (since M 2 may correspond 
to A, B, or C), M 3 with probability 0.375, and M 4 with probability 0.125. 

If Mi occurs, the trend line must be A and the forecast error could 
be reduced to zero by setting «i = % and /?i = %. Similarly, if M 4 occurs 
the trend line must be D; again «i = % and fa = Va are appropriate. 

If M 2 occurs, the trend line is either A, B, or C with equal probability. 
Since B is halfway between A and C in year 2, it follows that the mean 



28 THE BELL SYSTEM TECHNICAL JOURNAL, JANUARY 1 982 



square forecast error is minimized by forecasting the level of B in year 
2; this result is obtained with ai = % and /?i = Vs. Similarly, if M3 
occurs, the forecast with these gains is on C in year 2, and the forecast 
error is minimized. 

Thus, the optimal gains are a\ = %, /?i = Va and, after processing the 
second data point, the mean square forecast error is 0.5(0) 2 + 0.25(2) 2 
+ 0.25(— 2) 2 = 2, which is identical to the initial value. Thus, for this 
example there is no reduction in forecast error after processing the 
second data point. 

4. 1.4 G unknown 

The above examples assumed that G was known exactly. Suppose 
now, however, that the gains are chosen under the assumption that 
G = oo, i.e., ai = fii = 1, but in fact G = 0. In this case, which is 
illustrated in Fig. 10, the mean square forecast error in year 2 is 0.25(3) 2 
+ 0.25(1) 2 + 0.25(-l) 2 + 0.25(-3) 2 - 5, which is five times the initial 
value. 

This result, although exaggerated, is important because it shows 
that spa could actually perform worse than the conventional forecast- 
ing procedure. Indeed, this fact, combined with the practical impossi- 
bility of estimating G without error, is the major consideration in our 
decision to use the robust version of spa described in Sections 4.2.3 
and 4.3.2. 



A 




•AM, 



B 




A MEASURED LOAD 

»■ FORECASTED LOAD 

TRUE LOAD 



1 2 

YEAR 



Fig. 10 — Nonoptimal filter: assumed G = oo, actual G = 0. 
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4.2 Linear growth 

4.2.1 G known 

The results of this section apply when the true load has constant 
incremental growth and when the ratio G, eq. (21), is known exactly. 

Figure 11 displays the normalized rms 1-year forecast error as a 
function of the number of data points processed after initialization. 
The normalization is with respect to the rms forecast error of the 
conventional projection method. Since this method is used to initialize 
spa, the initial normalized rms error equals 1.0. Each curve corresponds 
to a different value of G and, accordingly, to a different set of gains 
ct n and 0„. For example, the gains corresponding to three different 
values of G are shown in Fig. 12. 

The results shown in Fig. 11 assume that the true load displays 
constant incremental growth; under this assumption, the forecast error 
approaches zero as n increases. In practice, however, unexpected 
changes in trend will occur and, consequently, as described in Section 
2.3, spa will occasionally be reinitialized. Accordingly, average forecast 
error over the interval between reinitializations is a more meaningful 
figure of merit. In the following paragraphs, we assume an average 
interval of five years. We emphasize, however, that this assumption 
will underestimate the benefits of spa under more stable conditions 
and overestimate the benefits under less stable conditions. 




Fig. 11— Optimal filter: forecast error vs. year. 
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Fig. 12 — Optimal gains vs. year. 

Thus, the curve labeled optimal in Fig. 13 displays the normalized 
5-year average forecast error — the average of the first five values in 
Fig. 11 — as a function of the parameter G. If, for example, G = G2 and 
we choose the gains accordingly, the average rms error would be about 
20 percent less than that for the conventional forecasting method. 
Alternatively, if G = G4, we would choose a different set of gains and 
the average reduction would be approximately 15 percent. Note as 
suggested by the examples discussed above, that the reduction in 
forecast error is a convex function of G. 

4.2.2 G unknown 

Suppose now that the actual value of G differs from the assumed 
value. For example, suppose we use the gains corresponding to G = G2, 
but in fact G = G 4 . In this case, the curve labeled G = G% in Fig. 13 
shows that the 5-year average rms error would be about 20 percent 
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Fig. 13— Average rms forecast error for several filters. 

larger than its initial value, in agreement with the example discussed 
in Section 4.1.4. 

The results of Fig. 13 can also be interpreted as follows: Suppose 
that G = G 4 is appropriate when trunk group loads are derived from 
trunk group data. Then G = G 2 = G 4 /3 would be appropriate when 
loads are derived from point-to-point (cmds) data since, from Fig. 2, 
the rms measurement error for cmds data is approximately three times 
that for trunk group data. Figure 13 then shows that an spa tuned for 
point-to-point data would perform poorly with trunk group data. 

The results shown in Fig. 13 might suggest that we should use 
different gains for different applications. We considered this approach 
but found it to have two practical problems. 

First, as shown in Fig. 2, measurement error depends not only upon 
the type of data but also upon load size. And, although we could allow 
the gains to be a function of load size and data type, multiple versions 
of spa would be relatively more difficult to implement and maintain. 

Second, and more important, in practice it will not be possible to 
determine the exact value of G. For, in general, actual loads will not 
display constant incremental growth, but may exhibit random fluctu- 
ations about a trend line. In this case, the model described in Section 
III is still appropriate if x n denotes the trend level, instead of the true 
load. However, a 2 now consists of two components: one because of the 
difference between the measured and true loads and the other because 
of the difference between the true load and the trend line. Although 
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Ref. 2 quantifies the first component, we have no way a priori to assess 
the contribution of the second component. Thus, our assumed value 
for a 2 may differ from the actual value. 

Similarly, since we cannot determine a priori the difference between 
our estimate of the aggregate and individual trunk group load growth 
rates, our estimate of o%, eq. (19), will, in general, differ from the actual 
value. 

Consequently, the assumed value of G will in general also differ from 
the actual value. Accordingly, to guard against the consequences of an 
error in our estimate of G, we would like to employ an algorithm which 
performs well under a variety of possible operating conditions. 



4.2.3 Robust gains for linear growth 



Fortunately, there exists a robust set of gains which provides the 
same average performance independent of the true value of G. That is, 
as shown by the curve labeled "robust" in Fig. 13, if we use the gains 
corresponding to G = G 3 (see Fig. 12) then the 5-year average rms 
error will be about 10 percent less than that of the conventional 
projection method — independent of the actual value of G. 

Of course, as with any robust technique, Fig. 13 shows that we pay 
a premium by receiving less than the theoretically optimal perform- 
ance for protection against the possibility of a performance substan- 
tially worse than that of the conventional projection method. However, 
we used the data from our field test 5 to estimate, a posteriori, the value 




Fig. 14 — Robust filter (variable gains): rms forecast error vs. year. 
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of G. Remarkably, the observed value was G » G 3 . Thus, at least for 
this one case, our robust filter was in fact nearly optimal. 

Although the robust gains yield an average performance which is 
independent of G, Fig. 14 shows that the actual performance in each 
year is not independent of G. For example, if G = G 5 the rms error of 
spa increases above that of the conventional method by about 10 
percent in year 2, but is decreased by about 35 percent in year 5. Thus, 
during the first couple of years after initialization, spa may perform 
slightly worse than the conventional method for some trunk groups. 
Thereafter, spa will perform better, provided that the trend remains 
constant. 

4.3 Robust, constant gains for SPA 
4.3. 1 Response to unexpected changes in trend 

With the "robust gains for linear growth," spa is robust to measure- 
ment and initial growth estimation errors. For practical applications, 
however, it is equally important that spa also be made robust to 
deviations from linear growth. 
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Fig. 15— Optimal gains with modelling error. 
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That is, the results of Section 4.2.3 assume that the trend line 
displays constant incremental growth. Under this assumption, since 
the forecast error approaches zero, the gains a n and /?„ approach zero 
as n increases. Accordingly, as discussed in Section 2.3, spa would 
eventually respond to unexpected changes in the level or slope of the 
trend only when their cumulative effect produced two consecutive 
outliers of the same sign. 

4.3.2 Constant gains 

At the expense of receiving less than the theoretically optimal 
performance in the ideal case where the trend remains constant, we 
can decrease spa's response time to unexpected changes in trend by 
not allowing the gains a„ and /?„ to approach zero. As discussed in Ref. 
7, one approach to selecting limiting values for a„ and /?„ is to add 
zero-mean, uncorrected random variables (w n and w n ) to the descrip- 
tion of the true load in eq. (9). That is, we assume that the covariance 
matrix Q of modelling errors is nonzero. These modelling error terms 
lead to gains a„ and /?„ which approach nonzero constants that depend 
upon the mean square value of these error terms. For example, Fig. 15 
shows the optimal gains corresponding to G = G3 for several values of 
E(wl) and E(wl). 

Although theoretically appealing, the above approach leaves the 
practical problem of determining appropriate values for E(Wn) and 
E(ibn). Values could be gleaned from a large amount of historical data, 
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Fig. 16— Robust filter (constant gains): rms forecast error vs. year. 
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which we do not have, but there is no guarantee that they would be 
appropriate for the future. 

Therefore, instead of pursuing the above approach, we simply re- 
placed the variable gains corresponding to G = G 3 in Fig. 12 by the 
average of the first few terms and made the following observations: 
First, as shown in Fig. 16, the performance in the ideal case where the 
trend remains constant is nearly identical for the first six years after 
initialization to that for the robust gains for linear growth. Thereafter, 
the performance is somewhat poorer, but we anticipate that only a 
small fraction of groups will remain on a constant linear trend beyond 
6 years. Second, as illustrated in Fig. 17, when the true load displays 
random deviations from a linear trend, i.e., when w n , w n are nonzero, 
the constant gains perform better than those designed for linear 
growth. Thus, although the performance is somewhat degraded in the 
ideal case, the use of constant gains provides protection against the 
very real possibility of unexpected changes in trend. For comparison, 
Fig. 17 also displays the performance with gains obtained by truncating 
the robust gains for linear growth in year 6. As indicated, performance 
is somewhat better in the ideal case, but worse when modelling errors 
are present. Since we do not have sufficient historical data to estimate 
the level of modelling errors, we are, therefore, recommending the use 
of the constant gains. 
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Fig. 17— Forecast error with modelling error. 
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Fig. 18 — Root-mean-square observed forecast error vs. load: outlier thresholds. 



V. OUTLIER THRESHOLDS 

As discussed in Section 2.3, Huber's studies 5 show that the outlier 
thresholds may be set anywhere in the range of one to two times the 
rms observed forecast error. In theory, the observed forecast error 
depends upon the number of data points processed after initialization; 
however, since the average rms forecast error decreases by only about 
10 percent and since performance is not sensitive to the precise setting 
of the thresholds, we will set them based upon the initial value of the 
rms observed forecast error; that is, by 

p 2 = E(x 1J0 - yi ) 2 . 



Since Xyo = (1 + g)yo and Xi = (1 + g)xo, it follows that 

p 2 = E[{y - x ) + ( v - x )g + x (g - g) - ( v, - x,)] 2 
= 2<r 2 (l +g + §*/2) + s 2 o 2 g , 
or, since |^| <c 1, 

p as xoOg + Zo . 



(22) 



(23) 



(24) 



To establish numerical values for p, we use the theoretical expres- 
sions for mean square measurement error given in Ref. 2: 



a 2 = 



20 



2x h 



+ V 



(25) 



where h is the mean call-holding-time in hours, p is the sampling rate 
(for trunk group data, p = 1.0; for cmds data, p = 0.05), and 



V d = max(0, 0.13*;? - 2x h) 



(26) 
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is the variance of the daily source loads. For first-routed loads, <J> = 1.5 
is usually appropriate. 

Figure 18 displays p as a function of the load x for both trunk group 
data (p = 1.0) and cmds data (p = 0.05). In each case, we assume 
o K = 0.06, which is the value observed during our field test, 6 and h = 



X A 2 . 



Based upon the results of Fig. 18 and Huber's studies, 5 it follows 
that outlier thresholds set at two times the value of p for trunk group 
data should provide adequate performance for both trunk group and 
cmds data. That is, this setting places the thresholds within the 
allowed range of about one to two times the actual value of p for both 
trunk group and cmds data. 

VI. CONCLUSIONS 

The spa employs a linear two-state Kalman filter, together with 
logic for detecting and responding to outlier measurements. The pa- 
rameters of spa have been selected to provide improved performance 
over the range of operating conditions, including the use of either 
trunk group or point-to-point traffic measurements. 

Under the assumption of linear growth for 5-year intervals, the 
average rms 1-year forecast of spa is about 10 percent less than that 
of the forecasting methods currently in use in the Bell System. More- 
over, the filter parameters and outlier procedures have been designed 
so that spa will respond to changes in trend. 

Field test results confirm the theoretical results presented here. 
Accordingly, specifications have been written for the inclusion of spa 
in the Bell System's standard mechanized trunk forecasting systems. 
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